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Anisotropic grid adaptation is examined by decomposing the steps of flow solution, ad- 
joint solution, error estimation, metric construction, and simplex grid adaptation. Multiple 
implementations of each of these steps are evaluated by comparison to each other and ex- 
pected analytic results when available. For example, grids are adapted to analytic metric 
fields and grid measures are computed to illustrate the properties of multiple independent 
implementations of grid adaptation mechanics. Different implementations of each step in 
the adaptation process can be evaluated in a system where the other components of the 
adaptive cycle are fixed. Detailed examination of these properties allows comparison of 
different methods to identify the current state of the art and where further development 
should be targeted. 


I. Introduction 

Output-based grid adaptation^ methods have been developed to control the discretization error of nu- 
merical simulation. Ficlkowski and Darmofal 1 provide a review of existing methods that encompass many 
different discretization schemes. Surprisingly, there are few examples where these methods have been di- 
rectly compared to each other on the same geometry. One example of this comparison is when anisotropic 
tetrahedral and Cartesian cut-cell methods were applied to predict near field sonic boom. 2 

Output-based grid adaptation is an iterative process that involves multiple steps, see Fig. 1. When the grid 
adaptation process is separated into these constituent parts, different methods for error estimation, metric 
creation, and grid adaptation mechanics can be compared at a finer grain level than simply monitoring the 
convergence of the output with grid optimization. This work will focus on adaptation of simplicies (triangles 
and tetrahedra) to an anisotropic metric. The specific details of the adaptive process have been examined 
in a few situations (e.g., Hessian recovery, 3 the formulation of the error estimate.) 4 

Often, output-based approaches are examined by comparing an adapted discrete solution to an analytic 
solution when available. To address more complex problems without a closed-form solution, an adapted 
solution can be compared to a solution computed on a much finer uniformly refined grid. While these 
evaluation methods are good global indicators of the correctness and efficiency of the entire output-based 
process, they are unable to expose weaknesses or inefficiencies in the individual steps of the process. 

This paper describes an approach targeted at improving the robustness and efficiency of the grid adap- 
tation process by decomposing the process and comparing different implementations of each step (i.e., each 
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^output-based, goal-based, and adjoint-based grid adaptation are common synonyms for each other. 
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Figure 1. Output-based grid adaptation process. 


box in Fig. 1) to understand the implications of implementation choices. Once the process is decomposed, 
each step can be verified by comparison to alternative implementations or an analytic result when available. 
Descriptions of each implementation and the results of this comparison are documented to ensure correct 
implementation and set the stage for further method development. This verification by comparison approach 
is also employed by the Turbulence Modeling Resource Website. 5 “What makes the current website unique is 
that it focuses on providing ready access to equations, grids, and flow solution details from previously verified 
codes as an aid to users who wish to verify their own implementations of models on relatively simple cases.” 6 
The goal of this work is to define a framework for rigorous examination of the entire output-based adapta- 
tion process and each of its constitutive parts that can guide the implementation and further development 
of solution adaptive methods. 


II. Definitions 


A metric tensor M. is commonly employed to define the desired multidimensional grid resolution. This 
symmetric positive definite matrix can be diagonalized: 


M = XAX t , 

where the eigenvectors X define an orthonormal basis with length specifications in this basis, 
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The metric can be interpreted as an ellipsoid with axes of lengths hi in directions of Xj. 7 

In a continuous setting, M is defined throughout the domain Q and the edge length I of an edge param- 
eterized between nodes A and B is given as, 


i M (A,B) = 


AB T M((1 - t)A + tB)AB d t, 


( 3 ) 


where the parameter t varies linearly between t = 0 at A and t = 1 at B. In this study, A4 is sampled 
at each node and interpolated to evaluate edge lengths (and shape measures used internally by the grid 
adaptation mechanics). The length of an edge in the metric £ is either computed with quadrature or 
geometric interpolation as given by Alauzet. 8 With a metric Xi defined at each node in the mesh, the length 
of an edge defined between nodes A and B is 

t M (A, B) = ,Jab t M(A)AB, (4) 


• In / 


where 


As r 


y/ AB T M(A)AB 
■sj AB T M{B)AB ' 

1, the A4(A) and A 4(B) become equivalent. For |r — 1| < 10 -12 , 

y/ AB T M(A)AB + \J AB T M(B)AB 

£m{ a iB) = , 
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to prevent the rlnr term from producing a divide-by-zero in finite precision arithmetic. The goal of adap- 
tation is to produce a grid with all edges of unit length 1=1. 

The statistical moments are computed in log space by 

1 E 

£= ^^ l0g 2^)’ ( 7 ) 

i-1 

where A is the length of edge i in AT There are E edges in the grid. The mean edge length reported in the 
table of statistics is 2 A The higher order moments are given by 

At = ^) X! l lo S2 (4) ~ • (8) 

^ i=l 

A series of adapted anisotropic grids do not have the a uniform scaling of element sizes, which are present 
in uniformly refined grids. For adapted grids, the complexity C can be used as a continuous counterpart of 
the classical parameter h from uniformly refined grid studies. 9,10 In the continuous setting, the complexity 
of a grid is defined as 

C= f v'det(JW) dll, (9) 

Jn 

for the domain O. In the discrete setting, M is sampled at each node i and 

N 

C = Y. ] \/det (Mi)Vi, (10) 

i-l 

where V) is the volume of the Voronoi dual surrounding each node. The complexity has a linear dependency 
with respect to the number of nodes and elements. It is used in error estimates to specify the desired accuracy 
via the size of the adapted grid. 

III. Output-Based Grid Adaption 

The output-based grid adaption cycle described in Fig. 1 is composed of various steps detailed in this 
section. Given an initial grid, the process starts by computing a flow solution. 

III. A. Flow and Adjoint Solution 

FUN3D 11,12 is used to compute the flow solution. An exact discrete adjoint 13 is computed with a solution 
scheme described by Nielsen et al. 14 that is wrapped with the Generalized Conjugate Residual (CGR) 15 
iterative scheme to accelerate convergence. Some of the viscous cases required the Hierarchical Adaptive 
Nonlinear Iterative Method (HANIM) to reach a converged steady-state residual on the high aspect ratio 
adapted grids. (HANIM was formerly referred to as HANIS by Pandya et al. 16 ) 

III.B. Error Estimation and Metric Construction 

Two error estimation, localization, and metric construction methods are compared in this study. These 
include the anisotropic optimal goal-based (Opt-Goal) metric of Loseille, Dervieux, and Alauzet 17 and the 
output-based anisotropic metric of Venditti. 18 Both of these methods require Hessian recovery 19 from the 
discrete solution. 

Hessian recovery is performed with the “double L 2 -projection” method of Alauzet and Loseille. 3 This 
method computes the gradient in each simplex based on the solution at the nodes. Then this reconstructed 
gradient is volume weighted to form an average reconstructed gradient from surrounding simplicies to each 
node. A second pass of this gradient reconstruction operator is used to recover second and mixed partial 
derivatives. Off-diagonal Hessian terms are averaged to yield a symmetric Hessian. 

Vallet et al. 20 provides analysis that shows that the reconstructed gradient is less accurate on the bound- 
aries than the interior. To prevent the inaccurate reconstructed boundary Hessian from negatively impacting 
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the adaptive process, the boundary values of the reconstructed Hessian are discarded and the boundary Hes- 
sian is averaged from interior nodes connected to the boundary nodes by edges. This averaging technique 
dramatically improved results and the convergence of the adaptive process. 

Details of the 3D implementation of the Venditti 18 metric are provided by Park . 21-23 Anisotropy of 
the metric is set by the Mach number Hessian, which is scaled so that its largest eigenvalue matches a 
spacing derived from the current spacing and an adjoint error estimate. This current spacing is estimated 
by computing the metric of each simplex 24 and averaging the implied metric to the nodes in log-Euclidean 
space. At the nodes, the largest eigenvalue is computed to estimate the local spacing of the grid. This 
process can produce noise in the current spacing estimate and is negatively impacted by sliver cells. So, 
two passes of isotropic gradation control are used to limit the local spacing to not less than 3/2 times its 
neighbor’s spacing. This prevents an isolated sliver element from indicating an inappropriately small local 
spacing. 

The Opt-Goal metric has the benefit of no direct dependency on the current grid, making it less sensitive 
to the size and shapes of the current grid elements. This metric depends only on the Hessian of the Euler 
fluxes weighted by the gradient of the adjoint state, 

•^Opt-Goal = (det |.ffgoal|) 2p+d |-Hgoal|i ( 11 ) 

with the norm order of p = 1, the dimensionality d, and the goal based Hessian, 

5 

tfgoai = H ([Aa^fx) + [Aj/]j(x) + [Az]j(x)) . (12) 

i= i 
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with W* denoting the j th component of the adjoint vector W* and H (TjiW : j)) the Hessian of the j th 
component of the vector Ti(W). 

The aspect ratio dR of the metric can be optionally controlled by replacing hi with hi and with /13 in 
Eq. ( 2 ): 

hi = max(/i 2 , dRdi) (14) 

/13 = max(/i 3 , jRhi). (15) 

Gradation control is implemented in the form of Eq. (9) from Alauzet 8 and applied to both solution adaptive 
metrics. A gradation 8 of (3 = 1.5 is used for all cases. The metric is rescaled to maintain the desired 
complexity after each of ten passes of gradation control. This is important for constructing a metric that is 
realizable by the grid adaption mechanics. 


III.C. Grid Adaptation Mechanics 

The mesh mechanics from four grid adaption tools are considered: BAMG 25 (Bidimensional Anisotropic Mesh 
Generator), refine, 26-2 ' Feflo.a , 28 and EPIC 29 (Edge Primitive Insertion and Collapse). Each of these 
codes attempts to produce a simplex (triangles in 2D and tetrahedra in 3D) mesh from a supplied anisotropic 
metric field. The focus of this study is to compare grid mechanic implementations at the elemental operator 
level. Variations in implementation near curved boundaries is avoided by only examining geometry models 
with linear surfaces. The impact of curved boundaries will be investigated in a future study. Results for 2D 
and 3D test cases with prescribed and computed sizing metrics are presented. 

The BAMG 25 (Bidimensional Anisotropic Mesh Generator) is used for 2D examples. The original version of 
refine, refine/one , 26 has been modified to target edges with a length in metric space outside of [0.3, 1.8]. 
The original intent 26 of the tool was to always produce edges that are equal or smaller than the metric 
specification. This study is focused on metric conformity or producing a unit mesh in the metric space, 
which requires a distribution of edge lengths centered about unity. The metric supplied to the tool is used 
as provided with no modifications (e.g., smoothing, limiting). A second version, refine/two, 2 ' is also 
evaluated, which only utilizes split and collapse operators . 29 The 2D algorithm has been extended to include 
the vertex smoothing of Alauzet . 30 The metric is used by refine/two as provided with no modifications. 
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Feflo.a belongs to the class of metric-based local remeshing software. It uses the classical operators 
of insertion, collapse, swap, and smoothing described by Frey and George. 28 Surface and volume meshes 
are adapted in a couple of ways to ensure a high-level of robustness. It is a based on a unique cavity- 
operator, where each of the standard operators (e.g., insertion, collapse, swap, smoothing ) are recast within 
an insertion procedure. The operation that results from this insertion depends only on the choice of the 
cavity. Cavity optimizations (e.g., reduction, enlargement) are performed to improve the quality of mesh 
after the insertion. Feflo.a handles 2D, surface, and 3D meshes. 

Grids from the EPIC adaptation code are presented for the 3D examples. The grid mechanics in EPIC have 
been extended to include element reconnection (swapping) and node smoothing operators. Element swapping 
is performed using the generalized swapping technique of Alauzet. 30 Node movement is performed using a 
local optimization technique to maximize a metric-based cell quality objective function. Edge lengths in 
EPIC are computed with numerical integration of Eq. (3) instead of assuming a linear variation of the metric 
M in log-Euclidean space. EPIC grids using three sets of mesh operators are presented: only insertion and 
collapse (EPIC-IC); insertion, collapse, and swap (EPIC-ICS); insertion, collapse, swap, and node movement 
(EPIC-ICSM). 

IV. Prescribed Metric Fields 

To examine the grid adaption step of Fig. 1 in isolation, 2D and 3D versions of an analytic metric field 
are used to evaluate the grid adaptation mechanics on a unit Cartesian domain in the following subsections. 

IV. A. Analytic 2D Metric Field 

The BAMG, refine/two, and Feflo.a grid adaptation processes are applied to an analytically defined (on 
the domain x = [0, 1] and 2 = [0, 1]) 2D metric field given by 

M = ( ) where 

V o h 2 ) 

The metric is provided at the nodes of an initial coarse grid. The grid is adapted to the metric on this 
grid and then the analytic metric is evaluated on the new grid. This process, where the analytic metric is 
evaluated at the nodes and the grid is adapted, is repeated nine times. By this point, the grid adaptation 
cycle has reached steady state. The final grids are shown in Fig. 2. Histograms of edge length i in A4 
evaluated on the final grids are provided for these methods in Fig. 3 and statistics on these final grid edge 
lengths are shown in Table 1. Distributions in Fig. 3 that appear more narrow and peaked produce lower 
values of the moments in Table 1. 

The minimum and maximum edge lengths are comparable for all methods, and bounded by the range of 
refine/two between 0.601 and 1.472. BAMG has a slightly smaller mesh with a vast majority of the edges near 
unit length in the metric. The methods produce grids with node counts of 0.597 to 0.723 of the complexity 
for this 2D metric. The theoretical minimum for the ratio of nodes to complexity is 

= 0.45345, (17) 

2 i/T2 V 

where 1/2 is the ratio of nodes to triangles in an infinite 2D grid with an average node degree of six and 
7r/-\/l2 is the density of triangles in an ideal Delaunay grid (also true for a constant metric field and ideal 
mapped Delaunay grid). 31 BAMG produced a grid with a node count that is closest to this theoretical limit. It 
also had the best statistics in terms of mean and higher order moments, which can be observed qualitatively 
in the Fig. 3 histograms. 
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Table 1. 2D linear metric properties. 


Mechanics 

Nodes 

N 

Complexity 

C 

Ratio 

N/C 

Min Edge 
min(t') 

Max Edge 
max(Q 

Mean 

Higher-order statistics 
M2 M3 M 4 

BAMG 

830 

1390.2 

0.597 

0.683 

1.436 

1.039 

0.0277 -0.0007 

0.0024 

refine/two 

1005 

1389.9 

0.723 

0.601 

1.472 

0.950 

0.0797 -0.0023 

0.0149 

Fef lo . a 

941 

1389.1 

0.713 

0.617 

1.413 

0.948 

0.0459 0.0026 

0.0055 



X 

(c) Feflo.a. 

Figure 2. Final adapted grids for the 2D linear metric. 
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(a) BAMG. 



(b) refine/two. 



(c) Feflo.a. 

Figure 3. Histogram of edge length in 2D linear metric. 
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IV. B. Analytic 3D Metric Field 


The Feflo.a, refine/one, refine/two, and multiple EPIC grid adaptation tools are applied to a 3D ana- 
lytically defined (on the domain x = [0, 1], y = [0, 1], and z = [0, 1]) metric field given by 

0 0 \ IX = o.i 

h~ 2 0 where } ly = 0.1 • (18) 

0 K 2 / [h z = 0.001 + 0.198 \z-0. 5| 

As in the 2D case, the metric is provided at the nodes of an initial coarse grid. The grid is adapted to the 
metric on this grid and the process is repeated until the process converges. The final grids are shown in 
Fig. 4. Histograms of edge length £ in A4 are provided for the two methods in Fig. 5. Statistics on these 
final grids are shown in Table 2. 

The ratio of the number of nodes to the complexity is near 2, which is expected by theory 8 9 and numerical 
experiments. 10 The Feflo.a grid is the smallest and has a large majority of near unity length edges in 
the metric space. Both versions of refine produce grids with more nodes than Feflo.a and a significant 
population of edges shorter than unity are shown in the histogram and minimum edge length in metric 
space. The ability of Feflo.a to prevent the creation of short edges or to remove existing short edges may 
be enhanced by the cavity operator. 32 The methods that include node movement (Feflo.a, refine/one, 
and EPIC-ICSM) tend to produce the histograms in Fig. 5 with a more peaked distribution. The methods 
without node movement (refine/two, EPIC-IC and EPIC-ICS) tend to produce a uniform distribution of 
edge lengths between their internal edge length limits. 

The addition of element swapping and node smoothing improves the statistics of the EPIC-IC algorithm. 
The long edges in the EPIC meshes can be attributed to the edge length computation method that EPIC uses 
compared to Feflo.a and refine. Statistics for the EPIC grid using the quadrature based edge lengths are 
shown in the last three lines in Table 2. When the edge lengths are evaluated using the quadrature based 
method consistent with the EPIC formulation, the maximum edge lengths become consistent with the other 
methods. 

Statistics for the EPIC edge lengths reported in Table 2 are computed with two methods. The first set 
are computed with Eq. (4), which is consistent with refine and Feflo.a. The second set with the footnote 
are computed using numerical integration of Eq. (3), which is consistent with EPIC’s internal calculations. 
The numerical integration approach is more accurate where the edge and metric principle directions are not 
aligned. The maximum and minimum edge lengths produced by EPIC are longer than refine. This reflects 
a design choice in the EPIC implementation that accepts a slightly longer edge if splitting that edge creates 
a very short edge. The addition of element swapping and node smoothing improves the statistics relative to 
EPIC-IC. 



Table 2. 3D linear metric properties. 


Mechanics 

Nodes 

N 

Complexity 

C 

Ratio 

N/C 

Min Edge 
min(t') 

Max Edge 
max(£) 

Mean 

2 M 

Higher- 

-order statistics 
M3 M 4 

Fef lo . a 

8933 

4668.6 

1.91 

0.449 

1.804 

1.011 

0.0597 

-0.0031 

0.0086 

ref ine/one 

9981 

4670.3 

2.14 

0.078 

1.799 

0.952 

0.1932 

-0.0506 

0.1341 

refine/two 

12868 

4673.8 

2.75 

0.065 

1.500 

0.977 

0.1423 

-0.0414 

0.1217 

EPIC-IC 

9524 

4671.4 

2.04 

0.149 

3.483 

1.024 

0.1227 

-0.0112 

0.0435 

EPIC-ICS 

9612 

4670.2 

2.06 

0.324 

3.051 

1.005 

0.1120 

-0.0054 

0.0314 

EPIC-ICSM 

9018 

4669.9 

1.93 

0.391 

2.439 

1.019 

0.0902 

-0.0094 

0.0215 

EPIC-IC * 

9524 

4671.4 

2.04 

0.149 

2.201 

1.022 

0.1209 

-0.0113 

0.0414 

EPIC-ICS * 

9612 

4670.2 

2.06 

0.324 

2.137 

1.004 

0.1113 

-0.0059 

0.0306 

EPIC-ICSM # 

9018 

4669.9 

1.93 

0.391 

1.907 

1.018 

0.0900 

-0.0095 

0.0214 

& These edge length statistics utilize numerical integration of Eq. (3), which is consistent with EPIC formulation. 

Other lines 


use Eq. (4). 
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(f) EPIC-ICSM. 


Figure 4. Final adapted grids for the 3D linear metric, 
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0.40 




(a) Feflo.a. (b) refine/one. 




(c) refine/two. (d) EPIC-IC. 




(e) EPIC-ICS. (f) EPIC-ICSM. 

Figure 5. Histogram of edge length in 3D linear metric. 
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V. Output-Based Adaptation 


Now that the adaptive mechanics has been evaluated with an analytic metric, the solution adaptive cases 
are examined. In this study, we consider only peace- wise linear geometries in 2D and 3D to avoid the surface 
approximation differences of these mesh generators. 

V.A. Diamond Airfoil in 2D Supersonic Flow 

A 2D diamond airfoil with a thickness of 0.07 is simulated in a supersonic inviscid flow with a Mach number 
of 2. The BAMG method is used to construct an initial coarse grid and a sequence of 20 grids with increasing 
complexity is computed with BAMG and refine/two. The complexity C 2 Diam and maximum aspect ratio 
A? 2 Diam are scheduled as 

C 2 oiam = min(1000c, 10000) (19) 

and 

A?2Diam = 10, (20) 

where c is the adaptive cycle index 

The final grids for the drag adaptation are shown in Fig. 6. The difference in drag of the final grid 
between the four methods is 0.00000174 (0.015 percent), see Table 3. The drag for this series of grids with 
increasing complexity are shown in Fig. 7. Drag computed on the INRIA Opt-Goal metric adapted grids 
approach the final drag more rapidly than drag computed on the Venditti metric adapted grids, and the 
BAMG method is able to produce grids that are closer to the requested complexity. 




0.6 
0.4 
0.2 
N 0 
- 0.2 
- 0.4 
- 0.6 


X 

(b) BAMG with Opt-Goal metric. 


x 

(a) BAMG with Venditti metric. 


0.6 
0.4 
0.2 
n o 
- 0.2 
- 0.4 
- 0.6 




(c) refine/two with Venditti metric. 


(d) refine/two with Opt-Goal metric. 


Figure 6. Final adapted grids for the 2D diamond airfoil. 
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Table 3. 2D diamond airfoil. 


Mechanics 

Error Estimation 

Drag 

C D 

Nodes 

N 

Complexity 

C 

Ratio 

N/C 

BAMG 

Venditti 

0.01133900 

10964 

10000 

1.10 

BAMG 

Opt-Goal 

0.01134039 

10518 

10000 

1.05 

refine/two 

Venditti 

0.01133865 

14136 

10000 

1.41 

refine/two 

Opt-Goal 

0.01134025 

14499 

10000 

1.45 



(a) Drag as function of degrees of freedom ( N ). 



(b) Drag as function of h = N 1 / 2 . 

Figure 7. Drag convergence for the 2D diamond airfoil. 
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V.B. Diamond Airfoil Extruded to 3D in Supersonic Flow 


A diamond airfoil with a span of one and a thickness of 0.07 is extruded to unit span in 3D supersonic 
inviscid flow with a Mach number of 2. The initial BAMG grid from the 2D case is extruded unit span to 
form two prisms from each triangle. This results in two symmetry planes and a center plane in the domain 
interior. This center plane enhances Hessian reconstruction accuracy, which degrades on boundaries. The 
extruded prisms are subdivided into the tetrahedra required by the simplex adaptation mechanics. The final 
grid sizes and drag values are shown in Table 4. The final symmetry plane grids for the drag adaptation 
are shown in Fig. 8. The drag for this series of grids with increasing complexity are shown in Fig. 9. The 
initial complexity is specified as 10,000 for the first ten adaptive cycles. This complexity is raised to 20,000, 
40,000, and 80,000 after every 10 subsequent adaptive cycles. Aspect ratio is not limited. The Venditti 
metric and the Opt-Goal metric are used for these cases. Feflo .a and refine/two are run sequentially, but 
refine/one is run on 8 cores with a domain decomposition scheme 26 and EPIC is run with 127 cores. 

As seen in Fig. 9, refine/one converged more slowly than the other methods to a specified grid size and 
drag. The Feflo. a method produces smaller grids for a given complexity and has a drag values closer to 
the 2D results. The EPIC-ICS method produces results similar to refine/two. There is a larger variation in 
drag (0.044 percent) between the adapted 3D grids than the 2D adaptive grids (0.015 percent). Finer grids 
could be attempted to verify that these methods are converging to the same drag value, because convergence 
to a common drag value is not obvious at this error level. Figure 9(b) indicates that the drag is increasing 
more rapidly for the Venditti error estimate adapted grids. As in the 3D analytic metric case, the ratio of 
the number of nodes to the complexity is near 2, which is expected by theory 9 and numerical experiments. 10 


Table 4. 3D extruded diamond airfoil. 


Mechanics 

Error Estimation 

Drag 

C D 

Nodes 

N 

Complexity 

C 

Ratio 

N/C 

Feflo . a 

Venditti 

0.01133797 

142026 

80000 

1.78 

refine/one 

Venditti 

0.01133649 

198679 

80000 

2.48 

refine/two 

Venditti 

0.01133657 

177023 

80000 

2.21 

EPIC-ICS 

Venditti 

0.01133698 

156656 

80000 

1.96 

EPIC-ICSM 

Venditti 

0.01133738 

154286 

80000 

1.93 

Feflo . a 

Opt-Goal 

0.01134097 

159740 

80000 

2.00 
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0.01135 r 


0.01134 - 


0.01133 - 


0.01132 - 


0.01131 - 


0.01130 



50000 


100000 

Nodes 


Venditti-Feflo.a 
Venditti-refine/one 
Venditti-refine/two 
Venditti-EPIC-ICS 
Venditti-EPIC-ICSM 
Opt-Goal-Feflo.a 
L 

150000 


200000 


(a) Drag as function of degrees of freedom (IV). 


0.01135 i- 


0.01134 - 


0.01133 


_L 


_L 


_L 


-e- 


-o- 


Venditti-Feflo.a 

Venditti-refine/one 

Venditti-refine/two 

Venditti-EPIC-ICS 

Venditti-EPIC-ICSM 

Opt-Goal-Feflo.a 



_L 


0.005 0.01 


0.015 


0.02 

h 


0.025 0.03 0.035 


0.04 


(b) Drag as function of h = N 1 / 3 . 

Figure 9. Drag convergence for the 3D extruded diamond airfoil. 
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V.C. Laminar Delta Wing (3D) 

The 3D Laminar Delta Wing case is described by Wang et al. 33 as part of its inclusion in all three of the 
International Workshops on High-Order CFD Methods to date. Adaptive results where originally presented 
by Leicht and Hartmann. 34 Drag and drag error are shown in Fig. 10. The reference drag used to compute 
drag error, (7j^ ference = 0.1656025, is the drag computed on a twice embedded final Feflo.a grid with 
10,853,343 nodes. The spacing h in the drag error plot is estimated using the suggestion of workshop 
problem description, h = iV -1 ' 3 . 

The same complexity schedule of the 3D extruded diamond airfoil is used in this example, which is 10,000, 
20,000, 40,000, and 80,000 for ten cycles each. The HANIM solution scheme is required to reach a converted 
steady-state residual on the two finest complexities. Feflo.a and refine/two are run sequentially, but 
refine/one is run on 8 cores with a domain decomposition scheme 26 and EPIC-ICS is run with 191 cores. 
Drag (Fig. 10) converges less rapidly in this case than the 3D Extruded Diamond Airfoil, which may be 
caused by the strong nonlinearities of this problem. The drag error is lower for a given number of degrees 
of freedom (DOF) than the errors of the higher-order methods summarized by Wang et al. 33 One possible 
explanation is the nominally second-order FUN3D method on adaptive grids may be resolving singularities 
that are not resolved by the workshop provided grids. 

Views of the final surface grid are provided for refine/one (Fig. 11), refine/two (Fig. 12), Feflo.a 
(Fig. 13), and EPIC-ICS (Fig. 14). These grids look very similar and clearly resolve the footprint of the 
leading-edge vortex on the upper surface of the delta wing. There is some slight clumping of the grid in the 
vicinity of the wing and symmetry plane intersection, which is suspected to be caused by noise from Hessian 
reconstruction. 

Overall, the final anisotropic adapted surface grids have a very similar appearance for all four methods. 
The drag error for this case is converging and is lower (for a given number of degrees of freedom) than the 
uniform grids of the High-Order CFD workshop. 
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(a) Drag as function of degrees of freedom. 



(b) Drag as function of h. 

Figure 10. Drag and drag error convergence for the Laminar Delta Wing. 
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(a) Upper wing surface intersection with symmetry plane. 
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(b) Upper wing surface. 

Figure 11. Final refine/one adapted grids for the Laminar Delta Wing. 
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(b) Upper wing surface. 

Figure 12. Final refine/two adapted grids for the Laminar Delta Wing. 
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(a) Upper wing surface intersection with symmetry plane. 
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(b) Upper wing surface. 

Figure 13. Final Feflo.a adapted grids for the Laminar Delta Wing. 
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(a) Upper wing surface intersection with symmetry plane. 
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(b) Upper wing surface. 

Figure 14. Final EPIC-ICS adapted grids for the Laminar Delta Wing. 


21 of 30 


American Institute of Aeronautics and Astronautics 



V.D. Turbulent Flat Plate (2D) 

This case is described on the Turbulence Modeling Resource 5 in the turbulence model numerical analysis 
section as the 2D Finite Flat Plate Validation Case. The geometry and boundary conditions are reproduced 
in Fig. 15. The Spalart-Allmaras (SA) 35 turbulence model is used for this case. This case has a length of 
symmetry boundary conditions before and after the no-slip plate to create more distance between these sim- 
ilarities and the inflow and outflow boundaries. This case has been examined with uniform grid refinement 36 
and adaptive grids. 3 ' 

The complexity Cpiate and maximum aspect ratio ^R p i a te are scheduled as 

Cpi a te = min(4000c, 80000) (21) 

and 

Opiate = 2°’ 5C , (22) 

where c is the adaptive cycle index. 

Drag convergence is shown in Fig. 16 for the first 23 grids. The series of adaptive grids appear to 
have a similar trajectory to CFL3D. This is suspected to be a coincidence, because other simulations (not 
shown) took a slightly different path. The Hessian recovery failed the on the 24th grid, computing not a 
number (NaN). The standard and 128-bit versions of BAMG is also applied (not shown), but failed on 8th 
grid adaptation when it was unable to find an enclosing triangle for inserting a new point. Skin friction is 
shown in Fig. 17 and profiles of velocity and eddy viscosity are shown in Fig. 18. The coarse isotropic initial 
grid and 23rd adapted grid are shown in Fig. 19. Zoom of grid near the end of the plate is shown in Fig. 20. 
The 3-axis as been scaled and A R is limited with Eq. (22), which makes the grid appear nearly isotropic in 
Fig. 20. This vertical scaling identifies triangles with lower aspect ratio (i.e., vertical stretching on the lower 
boundary) . 

The skin friction computed on the adapted grid is similar to the skin friction calculated on structured 
grids. The adapted grid skin friction is linear closer to the leading edge of the plate than the structured 
grid results, see Fig. 17(a) and (b). There is a small amount of noise in the adapted grid skin friction. 
The locations of the noise appear to correspond to the locations of the lower aspect ratio clumps on the 
lower boundary seen in Fig. 20. These may be a result of Hessian recovery issues. The w-velocity profile is 
indistinguishable from the structured grid results in Fig. 18(a). The largest deviation in the turbulent eddy 
viscosity profile is very close to the plate, see Fig. 18(b). The adaptive result shows no undershoot at the 
edge of the boundary layer. The adaptive scheme uses first-order convection for the turbulence model and 
the uniformly refined schemes use second-order convection. 

The small areas of grid clumping are suspected Hessian recovery issues on the lower boundary of Fig. 20. 
Without replacing Hessian recovery with averaging from the interior, inconsistent aspect ratios are requested, 
which produced grids with high amounts of noise in skin friction that stall the convergence of drag (not 
shown). This may be a contributing issue for the lower aspect ratio grid on the boundary exhibited by 
the equivalent FFMA method, which was reported by Yano and Darmofal. 38 Previous simulations did 
not converge the mean flow equations as well as the current simulation presented here. The remaining 
iterative error in these simulations is amplified by the Hessian recovery method and produces elements with 
inappropriate stretching and orientation in the interior (not shown). The “double L 2 -projection” Hessian 
recovery method produced less noise than applying an unweighted least squares gradient reconstruction 
method 11 twice (not shown). 
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Flat Plate Boundary Conditions, 
M=0.2, Re L =5 million (L=1), T ref =540 R 
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Figure 15. Turbulent flat plate geometry description. 



Figure 16. Turbulent flat plate drag convergence. 
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mu 




(b) Turbulent eddy viscosity profile at X = 1. 

Figure 18. Turbulent flat plate profiles. 


25 of 30 


American Institute of Aeronautics and Astronautics 



(a) Upper wing surface intersection with symmetry plane. 



(b) Upper wing surface. 

Figure 19. Turbulent flat plate initial and final grids. 
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Figure 20. Turbulent flat plate final grid (zoomed and scaled vertically). 


VI. Observations 

Having a metric field that is realizable (can be reasonably approximated by a discrete grid) is critical 
to robustness of grid adaptation mechanics. Having a metric field that describes a new mesh that better 
approximates the continuous solution is critical to robustness of the entire adaptive process. The adaptive 
process can be made more robust to poorly sized or shaped elements by eliminating the direct dependency 
on the current grid spacing or the metric of the current grid simplicies (i.e., Opt-Goal metric). 

Gradation control reduces discontinuous features and noise in the metric to improve the realizablity 
of the metric field. Hessian reconstruction, particularly on the boundary of the domain, is important for 
constructing a high quality metric field. To alleviate the known reduction in Hessian reconstruction accuracy 
on the boundary, the boundary Hessian reconstruction is replaced with averaging from the interior. This 
boundary reconstruction created a need for a center node in the extruded airfoil 3D case for averaging the 
Hessian to the symmetry planes. 

A strong iterative solver for the flow and adjoint solution is critical for the robustness of the complete 
adaptation cycle. The flow solution must be converged well or iterative error can be amplified by Hessian 
recovery. When this Hessian is used to construct a metric, elements with low aspect ratio and inappropriate 
orientation can be created for the turbulent flat plate case. 

The adaptive methods evaluated have ranges defined for acceptable long and short edges in the metric. 
These methods modify the grid in the neighborhood of these edges that are outside this acceptable range to 
produce a grid that conforms to the metric specification. The neighborhood of acceptable length edges are 
not modified. The choice of the acceptable length range in each scheme affects the statistics of edge length 
in the resultant adapted grid. 
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The implementation of the Opt-Goal metric and gradation limiting in FUN3D required detailed discus- 
sions with the authors of these methods. The availability of the test cases and results from this paper could 
make the implementation task easier and confirm correctness for these implementations. A database of grid 
adaptation examples would encourage the application of the Scientific Method 39 to speed the development 
of solution adaptive methods. Development of both turbulence modeling 6 and flight vehicle simulations 40 
has been significantly impacted by this approach, which justifies the expenditure of effort to build these 
databases. 


VII. Recommendations 

Further development of a verified, validated, fast, and robust flow and adjoint solvers is required for 
the use of high aspect ratio anisotropic grids to compute turbulent simulations. Further examination and 
development of Hessian recovery is needed (especially near boundaries and in the presence of iterative error) . 
This issue could be circumvented with a metric that does not require Hessian recovery. 4 

Gradation control implementations should be investigated by specifying a discontinuous metric and ex- 
amining the resulting grids after the application of gradation control. To extend the use of statistics to 
solution adaptive grids, the metric can be interpolated to the adapted grid and the edge length statistics can 
be computed. Statistics of the conformity 24 of the adapted grid to the metric field can also be computed 
with this interpolated metric. 

This paper has focused on simplex grid adaptation for planar geometries. There is an acute need for 
turbulent simulations of geometries with high curvature and complex topologies. New test cases should be 
constructed to address these needs (i.e., a more challenging analytic metric with high anisotropy intersecting 
a curved boundary). This paper utilized sequential execution with the exception of refine/one and EPIC. 
Evaluation of the scaling and correctness of parallel execution schemes for the grid mechanics are also required 
for large scale simulation. 


VIII. Conclusions 

An adaptive process is described and adaptive mechanics are evaluated by statistics of edge lengths in 
2D and 3D analytic metric fields. Euler, laminar, and turbulent solution adaptive cases are investigated 
with the adaptive mechanics. Important components of each step in the process are highlighted: flow solver 
robustness, Hessian recovery (particularly near boundaries), and gradation control. These concrete test cases, 
evaluated by each author, fostered detailed discussions about low-level implementation details between the 
authors. Hopefully, this paper will encourage follow-on papers, future organized paper sessions at conferences, 
a workshop, a website, or other dissemination of tests cases and results to spur future development of solution 
adaptive methods. 
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